Última actualización 23 de enero del 2021

Authores: Shalisko V., Castillo-Aja R., Santana E., Valdivia-Ornelas L.

Versión 6.1

Datos fuente

Datos abiertos SSA https://www.gob.mx/salud/documentos/datos-abiertos-152127

Fuentes auxiliares

Diccionarios de datos SSA https://www.gob.mx/salud/documentos/datos-abiertos-152127 y Catálogo Único de Claves de Áreas Geoestadísticas Estatales, Municipales y Localidades de INEGI https://www.inegi.org.mx/app/ageeml/

## Warning: package 'zoo' was built under R version 4.0.3
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Loading required package: sp
## rgdal: version: 1.5-15, (SVN revision 1045)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
## Path to GDAL shared files: C:/Users/vshal/Documents/R/win-library/4.0/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: C:/Users/vshal/Documents/R/win-library/4.0/rgdal/proj
## Linking to sp version:1.4-2
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.

ANALISIS MÉXICO

Pruebas por COVID-19 por fecha de inicio de síntomas

ruta_estados <- '../datos_shp/Estados.shp'
estados <- rgdal::readOGR(ruta_estados)
## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS =
## dumpSRS, : Discarded datum International_Terrestrial_Reference_Frame_1992 in CRS
## definition: +proj=lcc +lat_0=12 +lon_0=-102 +lat_1=17.5 +lat_2=29.5 +x_0=2500000
## +y_0=0 +ellps=GRS80 +units=m +no_defs
## OGR data source with driver: ESRI Shapefile 
## Source: "D:\GD\Projects_actual\geo_COVID-19\datos_shp\Estados.shp", layer: "Estados"
## with 32 features
## It has 2 fields
ruta_municipios <- '../datos_shp/Municipios.shp'
municipios <- rgdal::readOGR(ruta_municipios)
## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI, dumpSRS =
## dumpSRS, : Discarded datum International_Terrestrial_Reference_Frame_1992 in CRS
## definition: +proj=lcc +lat_0=12 +lon_0=-102 +lat_1=17.5 +lat_2=29.5 +x_0=2500000
## +y_0=0 +ellps=GRS80 +units=m +no_defs
## OGR data source with driver: ESRI Shapefile 
## Source: "D:\GD\Projects_actual\geo_COVID-19\datos_shp\Municipios.shp", layer: "Municipios"
## with 2456 features
## It has 4 fields
municipios@data$Clave_Mun_Ent_Texto <- as.numeric(as.character(municipios@data$CVE_MUN)) + 
                             1000 * as.numeric(as.character(municipios@data$CVE_ENT))
municipios@data$Clave_Mun_Ent_Texto <- sprintf("%05d", municipios@data$Clave_Mun_Ent_Texto)
#str(municipios@data)
metadatos_m <- datos_municipio[,1:12]
#str(metadatos_m)
dia_c_m <- datos_municipio[,c(casos_fechas)]
dia_n_m <- datos_municipio[,c(neg_fechas)]
dia_antes_c_m <- cbind(datos_municipio[,c(casos_fechas[-1])],data.frame(inicio = rep(0,dim(datos_municipio)[1])))
dia_antes_n_m <- cbind(datos_municipio[,c(neg_fechas[-1])],data.frame(inicio = rep(0,dim(datos_municipio)[1])))
incremento_m <- (dia_c_m + dia_n_m) - (dia_antes_c_m + dia_antes_n_m)
incremento_c_m <- dia_c_m  - dia_antes_c_m
names(incremento_m) <- pruebas_fechas
names(incremento_c_m) <- pruebas_fechas
#str(incremento_m)
incremento_m_a100k <- round(100000 * incremento_m / datos_municipio$Pob_Municipio, 2)
incremento_m[is.na(incremento_m)] <- 0 
incremento_c_m[is.na(incremento_c_m)] <- 0 
incremento_m_a100k[is.na(incremento_m_a100k)] <- 0

#dim(incremento_m_a100k)
#head(incremento_m_a100k)
semanas_lag <- c(seq(from = 0, to = n, by = 7))
pruebas_breaks <- c(0,10,20,30,40,50,75,100,150,1000000)
pruebas_breaks_labels <- c("0-10","10-20","20-30","30-40","40-50","50-75","75-100","100-150",">150")
pruebas_raw_breaks <- c(0,5,10,20,50,100,200,500,1000,1000000)
pruebas_raw_breaks_labels <- c("1-5","6-10","11-20","21-50","51-100","101-200","201-500","501-1000",">1000")
pruebas_indice_breaks <- c(-1,-0.75,-0.5,-0.25,-0.1,0.1,0.25,0.5,0.75,1)
pruebas_indice_breaks_labels <- c("menos que esperado (-1)","-0.75","-0.5","-0.25","esperado (0)","0.25","0.5","0.75","mas que esperado (1)")
pruebas_positividad_breaks <- c(0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,1)
pruebas_positividad_breaks_labels <- c("0-10%","10-20%","20-30%","30-40%","40-50%","50-60%","60-70%","70-80%",">80%")

#for (i in c(3,10,30,46)) {
for (i in 1:46) {
  selection_cols <- paste0("pruebas",fecha_cadena(as.Date(fecha_inicio) + seq(semanas_lag[i]+1, semanas_lag[i]+8)))  
  #print(selection_cols)
  suma_pruebas_semanal <- apply(incremento_m[,selection_cols], 1, sum)
  suma_pruebas_c_semanal <- apply(incremento_c_m[,selection_cols], 1, sum)
  
  total_pruebas_semanal <- sum(suma_pruebas_semanal, na.rm = TRUE)
  total_pruebas_c_semanal <- sum(suma_pruebas_c_semanal, na.rm = TRUE)
  total_poblacion <- sum(datos_municipio$Pob_Municipio, na.rm = TRUE)
  
  pruebas_positividad_semanal <- suma_pruebas_c_semanal / suma_pruebas_semanal
  
  pruebas_esperados_semanal <- total_pruebas_semanal * datos_municipio$Pob_Municipio / total_poblacion
  ## indice (similar a NDVI)
  pruebas_indice_semanal <- (suma_pruebas_semanal - pruebas_esperados_semanal) / (suma_pruebas_semanal + pruebas_esperados_semanal)

  suma_a100k_semanal <- apply(incremento_m_a100k[,selection_cols], 1, sum)
  datos_semanales <- data.frame(
             Clave_Mun_Ent_Texto = as.character(metadatos_m[,"Clave_Mun_Ent_Texto"]), 
             Municipio = metadatos_m[,"Nom_Mun"],
             Latitud = metadatos_m[,"Lat_Decimal"],
             Longitud = metadatos_m[,"Lon_Decimal"],
             pruebas_100k = suma_a100k_semanal,
             pruebas_raw = suma_pruebas_semanal,
             pruebas_esperado = pruebas_esperados_semanal,
             pruebas_indice = pruebas_indice_semanal,
             pruebas_positividad = pruebas_positividad_semanal
             )
  
  #print(str(datos_semanales))
  #print(head(datos_semanales, n = 30L))
  mis_municipios <- municipios
  mis_municipios <- merge(municipios, datos_semanales, by = "Clave_Mun_Ent_Texto", all.x = TRUE)
  mis_municipios@data$Class_100k <- as.numeric(cut(mis_municipios@data$pruebas_100k, 
                                                   breaks = pruebas_breaks))
  mis_municipios@data$Color_100k <- brewer.pal(n = 9, name = 'YlGn')[mis_municipios@data$Class_100k]
  mis_municipios@data$Class_raw <- as.numeric(cut(mis_municipios@data$pruebas_raw, 
                                                   breaks = pruebas_raw_breaks))  
  mis_municipios@data$Color_raw <- brewer.pal(n = 9, name = 'PuBuGn')[mis_municipios@data$Class_raw]  
  mis_municipios@data$Class_indice <- as.numeric(cut(mis_municipios@data$pruebas_indice, 
                                                   breaks = pruebas_indice_breaks))  
  mis_municipios@data$Color_indice <- brewer.pal(n = 9, name = 'PRGn')[mis_municipios@data$Class_indice]  
  mis_municipios@data$Class_positividad <- as.numeric(cut(mis_municipios@data$pruebas_positividad, 
                                                   breaks = pruebas_positividad_breaks))  
  mis_municipios@data$Color_positividad <- brewer.pal(n = 9, name = 'Reds')[mis_municipios@data$Class_positividad] 
  
  #print(mis_municipios@data$pruebas_positividad)
  #print(head(mis_municipios@data, n = 25L))
  #str(mis_municipios@data)
  
  #str(mis_municipios)
  
  ## grafica pruebas por 100 mil
  plot(mis_municipios,
          col = mis_municipios$Color_100k,
          axes = TRUE, 
          border = NA, 
          #col = "transparent",
          main = paste0("Pruebas semanales por 100 mil habitantes ",
                        as.Date(fecha_inicio) + semanas_lag[i] + 1," - ",
                         as.Date(fecha_inicio) + semanas_lag[i] + 8))
  plot(estados, add = TRUE)
  legend("topright", 
         c("Número de pruebas por COVID-19 con base en los reportes SSA de México"))  
  
  legend("bottomleft", 
         as.character(pruebas_breaks_labels),
         fill = brewer.pal(n = 9, name = 'YlGn')[1:9],
         title = "Por 100 mil habitantes"
         )

  ## grafica Pruebas
  plot(mis_municipios,
          col = mis_municipios$Color_raw,
          axes = TRUE, 
          border = NA, 
          #col = "transparent",
          main = paste0("Pruebas semanales ",
                        as.Date(fecha_inicio) + semanas_lag[i] + 1," - ",
                         as.Date(fecha_inicio) + semanas_lag[i] + 8))
  plot(estados, add = TRUE)
  legend("topright", 
         c("Número de pruebas por COVID-19 con base en los reportes SSA de México"))  
  legend("bottomleft", 
         as.character(pruebas_raw_breaks_labels),
         fill = brewer.pal(n = 9, name = 'PuBuGn')[1:9],
         title = "Pruebas semanales"
         )
  
  # ## grafica defuunciones indice
  # plot(mis_municipios,
  #         col = mis_municipios$Color_indice,
  #         axes = TRUE, 
  #         border = NA, 
  #         #col = "transparent",
  #         main = paste0("Índice de desviación del número de pruebas semanales esperado ",
  #                       as.Date(fecha_inicio) + semanas_lag[i] + 1," - ",
  #                        as.Date(fecha_inicio) + semanas_lag[i] + 8))
  # plot(estados, add = TRUE)
  # legend("topright", 
  #        c("Áreas en blanco reprezentan zonas sin pruebas realizados en el periodo",
  #          "Número de pruebas COVID-19 con base en los reportes SSA de México"))
  # legend("bottomleft", 
  #        as.character(pruebas_indice_breaks_labels),
  #        fill = brewer.pal(n = 9, name = 'PRGn')[1:9],
  #        title = "Indice de pruebas"
  #        )  
  # ## grafica defuunciones indice
  plot(mis_municipios,
          col = mis_municipios$Color_positividad,
          axes = TRUE,
          border = NA,
          #col = "transparent",
          main = paste0("Positividad en pruebas de COVID-19 ",
                        as.Date(fecha_inicio) + semanas_lag[i] + 1," - ",
                         as.Date(fecha_inicio) + semanas_lag[i] + 8))
  plot(estados, add = TRUE)
  legend("topright",
         c("Áreas en blanco reprezentan zonas sin pruebas realizados en el periodo",
           "Número de pruebas COVID-19 con base en los reportes SSA de México"))
  legend("bottomleft",
         as.character(pruebas_positividad_breaks_labels),
         fill = brewer.pal(n = 9, name = 'Reds')[1:9],
         title = "Positividad"
         )
  
  ## verification plot wit coordinates only
  #plot(x = mis_municipios@data$Longitud,
  #     y = mis_municipios@data$Latitud,
  #     col = mis_municipios@data$Color_raw,
  #     type = "p", pch = 19)
  
}

knitr::asis_output(htmltools::htmlPreserve(cc_html))
Creative Commons License.